ON-LINE TRACKING OF A SMOOTH REGRESSION 

FUNCTION 
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Abstract. We construct an on-line estimator with equidistant design 
for tracking a smooth function from Stone-Ibragimov-Khasminskii class. 
This estimator has the optimal convergence rate of risk to zero in sample 
size. The procedure for setting coefficients of the estimator is controlled 
by a single parameter and has a simple numerical solution. The off-line 
version of this estimator allows to eliminate a boundary layer. Simula- 
tion results are given. 



1. Introduction. 

In this paper, we consider a tracking problem for smooth function / = 
f(t)i < t < T, under observation 

X in = f{t in ) + ad, (1.1) 

for ti n = ^, i = 0, . . . , n (n is large), where is a sequence of i.i.d. 
random variables with E£i = 0, E£f = 1, and a 2 is a positive constant. 
Without additional assumptions on the function / it is difficult to create an 
estimator even for large re. The filtering approach, see Bar-Shalom and Li 
|J, proposes an estimator in the form of Kalman filter corresponding to a 
stochastic model for /, e.g. / is differentiable k times, and fc-th derivatives 
of / is simulated by a white noise with a certain intensity. Since / is deter- 
ministic function, the non-trivial part of such approach is a choice of filter 
parameters and asymptotic analysis of estimation risk in re — » oo. On the 
other hand, nonparametric statistic approach to the regression estimation 
of a function / assumes that / belong to some limited class. We take / 



from the class L) (introduced by Stone, [1C], [11| and Ibragimov and 
Khasminskii, ||, of k times continuously differentiable functions with 
Holder continuous last derivative (here /^°^ = /, L and a are the same for 
any function from the class): 

obeys k derivatives, /(°), . . . , f^; 
E(J3,L) = {f: |/W(t 2 )-/W(ti)|<L|t 2 -ti|« Vti,t2, a£(0,l]; 
(5 = k + a 

It is known |l(J, [11], [[|, [|| that there are kernel type estimators fn\t) of 
f^\t), j = 0, 1, . . . , k such that for a wide class of loss functions £,(*) and 
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p > 1 (C is positive constant): 

sup S£(n^||#-/W|| ip ) <C, j = 0,1,..., A; (1.2) 

/€S(/3,L) v y 

and no estimator provides a better rate of convergence to zero in n — > oo 
uniformly in £(/?, L). The same rate in n is valid under fixed value t for 

the estimation risk E(fn\t) — f^(t)) 2 , j = 0, 1 . . . , k. This rate cannot 
be exceeded uniformly on any nonempty open set from (0, T). Parallel to 
kernel type estimators (see, e.g. [§}]), ||), Khasminskii and Liptser @ 
proposed an on-line estimator (hereafter for brevity ij n is replaced by ij and 

# by fU)): 

n n 2/3+ 1 

j = 0,1,... ,k-l (1.3) 

f$(U) = f^Ku-i) + -$=s>{Xi - 7%-i)), 

subject to the initial conditions /(°)(0), fW (0), . . . , /W(0). The initial con- 
ditions are chosen as arbitrary bounded constants independent of n. The 
parameters qo, . . . ,q k are specifically chosen. The vector q with these entries 
is called the filter gain. 

The rigorous result given in Q is formulated as: 

Let filter gain q be chosen such that all roots of characteristic 
polynomial 

p k (u, q) = u k+l + q u k + qi u k ~ l + ... + q k _ lU + q k (1.4) 

are different and have negative real parts. Let the observation 
model defined in di.ij), f 6 S(/3, L) and a 2 > 0. Then, for the 

estimator given in ftl.Sj) there exist positive constants c(q), C(q) 

i 

(independent of n) such that for any ti > c(q)n 2/3+1 logn the 
normalized in n risk obeys 

k 

K^oo sup y j E{f^(t i )-7^{t i )fn^ <C(q). (1.5) 

2(/9-j) 

T/ie rates n 23+1 , j = 0, 1, . . . , k cannot be improved. The 



boundary layer c(q)n 2/3+1 logn, where (|1.5| ) might fail, is in- 
evitable. 

i 

Remark 1. T/ie Ze/fc side boundary layer c(q)n 23+1 logn is due to on-line 

limitations of the above tracking system. One can readily suggest an off-line 

modification with the same recursion in the backward time subject to some 

boundary conditions independent of observation X{ 's. This modification pos- 

i 

sesses the right side boundary layer [T — c(q)n 2/3+1 logn, T] and accuracy 
(\1.3j) on [0, T — c(q)n 2 i 3+1 logn]. So, some combination of the forward and 
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backward time tracking algorithms allows accuracy on [0, T] . For in- 

stance, a combination of forward time tracking on the interval and 
backward time tracking on [0, ^] can be used. 

In this paper, we deal with the estimator given in fll.3j ) and restrict our- 
selves by considering / from the class Y,(L,k + 1), i.e. the class of fe-times 
differentiable functions / with Lipschitz continuous f^ k \t). 

A suitable choice of filtering gain q should satisfy multiple requirements 
regarding the cost function C(q) and parameter c(q), involved in the descrip- 
tion of boundary layer. Moreover, a correct choice of q should guarantee that 
the roots of characteristic polynomial p k (u, q) are different and have negative 
real parts. These requirements might contradict each other. To avoid con- 
tradictions, we use the fact that estimator ( |1.3[ ) has a structure of Kalman 
filter. We build a Kalman filter according to Bar-Shalom and Li JlJ], so that 
f( k \t) is generated by a white noise with intensity 7. For each 7, we choose 
the Kalman gain q(7) and use it for minimization of C(q(7)) in 7. So, the 
minimization problem of the cost function is controlled by single parameter 
and allows to establish a reasonable relationship between C(q) and c(q). 
Moreover, this type of minimization automatically guarantees negative real 
parts of the roots for characteristic polynomial p k (u, q). For k < 4 the roots 
of p k (u, q) are different and numerical verification of the same fact for k > 4 
is available. 

2. The filter gain choice 



2.1. Preliminaries. Henceforth (3 = k + 1. For notational convenience we 
describe our problem in matrix notation. Introduce the following matrices: 

/o 1 ... o\ 
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\0 



1 



(fc+l)x(fc+i) 

M 



°)lxflH-lV b 

W (fe+l)xl 

Notice that by Lemma 3.1 in [2[ the roots of p k (u, q) and eigenvalues of 
(a — qA) coincide. In accordance with this remark, while eigenvalues of 
(a — qA) have negative real parts we may describe the cost function C(q) in 
terms of the bias M(q) and variance -P(q) for tracking errors (see, ( |3.11| )) 
(hereafter * is the transposition symbol): 

C(q) = trace (p(q) + M(q)M*(q] 

where M(q) = L(a — qA)~ 1 b and the matrix P(q) solves the Lyapunov 
equation (a - qA)P(q) +P(q) (a - qA)* + <r 2 qq* = 0. In Section we select 
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the filter gain q(7) from one-parameter family 

Q( 7 )A*~ 



{7 > le > : q(7) 



a 2 



where (5(7) is a positive definite matrix given by the algebraic Riccati equa- 
tion aQ(i) + Q{i)a* + 7 2 bb* - Q(7)A J 2 AQ(7) = 0. Finally we choose 

7 = argminC(q(7)) 

and the filter gain q(7°). A relevant choice of -y e allows to have an acceptable 
value of the constant c(q(7°)). 

2.2. Explicit formulae. In Section ||, we show that the cost function is 
expressed as: 

C« = a> ( jT ^^-^ + (§)'[(i)\g(a 

Furthermore, we give the explicit structure of ^(7) as a function of the 
control parameter (7 / 'a) l ^ k+l \ Namely 

/ry\ 1/k+l 

9o(7) = ^00 {-) 

/-y\ 2/k+l 

Qkh) = u ok[~), 

where Uij, i, j = 0, 1, . . . , k are entries of the matrix U being solution of the 
algebraic Riccati equation all + Ua* + bb* - UA*AU = 0. 



2.3. Example 1. Here, we consider the tracking problem for Lipschitz con- 

T 2 
' 00 



tinuous function /. Since k = 0, we have a = 0, A = 1, b = 1, 1 — C/q = 



and so, go = ^ Therefore, 



/ ^ °~7 l2 ° 2 

c(q(7)) = + 



2 7 2 • 

2/3 



With 7e < (2L)2/ 3 cr 1 / 3 , we have 7 = (2L) 2 / 3 a 1 / 3 and q(7°) = ^ 
The following estimator is constructed (here f(ti) '■= f^ n (ti)) 

f(U) = f(U^) + (—) (X t - f{U-i)). 



Remark 2. Notice that the direct minimization of C(q) = -P(q) + M 2 (q) 
with respect to qo provides the optimal qo = 9(7°) • -For k > 1, i/iis coinci- 
dence is not guaranteed. Under the direct minimization of the cost function 
C(q) with respect to q the eigenvalues of (a — (]*A) might have nonnegative 
real parts. 
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Figure 1 . Forward and backward time tracking with n = 2000 



2.4. Example 2. Let us consider a numerical solution for k = 2, i.e. / is 
twice differentiable function. Its second derivative is Lipschitz continuous 
with constant L = 100. For a = 0.25, we find (see Figure |) 7° = 24.533 
and 7°/cj = 98.132. According to Table 1, U 00 = 2, U 01 = 2, U 02 = 1. 
Hence qo = 9.225, q\ = 42.550, and q 2 = 98.132. So, the following estimator 
is constructed 

n n D /' 
- — - - — - i - — - 49 - 

7%) = mu-i) + - 7^(*i-i))- 

The combination of forward and backward tracking practically allows to 
eliminate the boundary layer (see Figure [l]) . 



3. Error analysis 



In this section, we present a vector of normalized tracking errors and 
derive the expression for C(q). 
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3.1. Notation. For notational convenience, set 



// (( %)\ 
Kf {k Ku)J 



and introduce a diagonal matrix C n and vector q n : 



g g_i g_fc 
C n = diag(n 2 / 3 + 1 ,n 2 / 3 + 1 ,... ,n 2 ' 3 + 1 ) and q n 



/ — 2fi_\ 
' g « 2/3+1 » 

_ 2/3-1 

qin 20+1 



2/3-fc 

\q k n ^ 



We use a vector-matrix form of estimator ( |1.3| ) 

n v 

and obvious identity 



(3.1) 



F(ti) = + -aF(t 8 _!) + b(/W(ti) - 



(3.2) 



3.2. Normalized errors. Denote <5j = -F(ii) — F{ti) and introduce normal- 
ized error Aj = C n <5j. Recursions ( |3.1[ ) and ( |3.2| ) provide 

<5j = <5j_i + -aSi-i + q n cr£i - q n -4<5j-i + (/(**) - /(*i_i))q n 
n 

-b(/ (fe) (ti)-/ (fc) (ti-i)). 
Multiplying both sides of this equation from the left by C n we find 

Aj = Aj_i -| — C n a5i-i + C„q n cj^ - C n q n A5,;_i 
n 

+ (/&) - /(ti_i))C n q n - 0»b - / (fc) (t,-i)). (3.3) 

A special structure (see 0) of the objects involved in (|3.3[) 



C n q 
1 



n 2 ^+!q 



2/3 



C n a5i-i = n 2 /3+iaA 



i-l 



2/3 



C n q n AS i - 1 =n w+iqAA^ 
C n b = nWif) 
allows to simplify fl3.3| ) significantly: 

2/3 /3 

A* = Ai_i + n 2 "+ 1 (a - q A) Aj_ x + n */>+i qcx£j 
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With D n = I + n 2 P +1 (a — qA) we rewrite the recursion for Aj's to 
At = D n Ai-i +n~W+-i<\cr£i 

+ n-W^q(f( tl ) - f(ti-x)) ~ nW^b{f^(t t ) - fW{U-l)) • (3-4) 

In Proposition 4.1 in Q, it is shown that for n large enough the magnitudes 
of all eigenvalues of D n are strictly less than 1. We shall use this property 
for asymptotic analysis and continuous time approximation. 

3.3. Normalized bias and variance. Denote the bias and variance of the 
normalized error A,;: 

Mf = EAi and ff = £(A 4 - Mf)(A; - Mf)*. 



Taking the expectation from both side of (3.4) we find 

= DnM^ (3.5) 



From Q and §), we get (A, - M") = L> n (A;_i - AfJL x ) + n 23+1 qCT £., 
so that Pf = E(Ai - Mf )(A, - Mf 1 )* is defined by the recursion 

PP = D n Pp^D* n + n~^a 2 qq*. 
Since - f^(U-i)\ < it is natural to choose / with 

n n 

and substitute supjgvv^) M"(M™)* by M?{Mf)*. Henceforth, is de- 
fined by recursion (|3.5|) with such /, i.e. 

~ gg , v 2/3 2/3 

Mf = M^_ 1 +n (a - qAjiW-Li +n g n - n ^Lb, 

where g n = jt^t^h ■ Thus, trace (P™ + M"(M J n )*) determines the normal- 
ized mean square tracking error 

C(q) =lmT ?woo trace + i\^ n (Mf )*) , > c(q)n"^ logn . ( 3 . 6 ) 



3.4. Continuous time approximation. To find "linin^oo" in ( |3.6D , we 
give a continuous time approximation of (M",P™). To this end, let us 
introduce the time stretching 

_1 2/3 

(*i — = n (si-*i_i)=n 

_ 1 

with to = s o = 0. The boundary layer [0, c(q)n 2 ' 3+1 logn] is transformed to 

[0,c(q)logn] and the interval [0,T] to [0,Tn*P+ 1 ]. Let us define M™ = Mf 
and Pg. = P™, i = 0, 1, . . . , and for s G Sj) 

= M«_ 1+ [ S f(a-qA)M2 i _ 1 )ds' + (e n -bL)(s-s^ 1 ) 

•ISi-X 

+n~ W (a - qA) (a - qA)* + a 2 qq*)ds'. (3.7) 
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. 1 , 

We shall consider these recursions for Sj from (c(q)logn, n 2 ' 3 + 1 Tj, where 

(M™,P™) have entries bounded in n (see, ||). Taking into account that 
recursions (|3.71 ) are homogeneous in s, let us replace {s and Sj} by {u = 

s — c(q)logra and ui = si — c(q)logn}. Then, the entries of (M",P™) are 

i 

bounded in n for < u < n 2 i 3+1 T — c(q)logn. Therefore without loss of 
generality we may consider (3/7) with initial conditions bounded in n: 

C(q) = lim Urn ™ trace (p% + M£(M£)*Y (3.8) 

u— >oo v v ' ' 



To determine the right hand side of (|3.8| ), we apply the Arzela-Ascoli theo- 
rem. For any T > 0, the functions (M",P™)o< u <t are uniformly bounded 
and equicontinuous. So, by the Arzela-Ascoli theorem, any converging sub- 
sequence (M™ , P™ ) obeys the limit (M^, P„) in the local uniform topology: 

lim ^ w min(l, sup (||M^' - M' u \\ + ||P U "' - P>J)) = 0, 

JV=1 0<u<N 7 

where 

M^ = M^ + j^ ((a-qA)M^ + bI?)dv 

Since the eigenvalues of a — qA have negative real parts, the limits M(q) := 
lim^oo M' u and P(q) := liniu^oo P' u exist and are defined as: 

M(q) = -L(a-qA)~ 1 b (3.9) 
(o - qA)P(q) + P(q)(a - qA)* + a 2 qq* = 0, (3.10) 

that is M(q) and P(q) are independent of {n 1 } and so 

C(q) = trace (P(q) + M(q)M*(q)) . (3.11) 

4. Minimization of the cost function in one parameter class 

4.1. Motivation. For large values of A: a direct minimization of C(q) from 
( p. 11 ) would be a difficult problem. Moreover, 

q° = argminC(q) 
q 

could not a priori guarantee negative real parts of eigenvalues for (a — q°A). 
To avoid implementation of a conditional minimization procedure, we pro- 
pose to choose q from some limited class given below. 

4.2. Adaptation to Kalman filter design. Our estimator has a structure 
of Kalman filter in the discrete time. We assume that P(0) is a random 
vector, P n (0) = EF(0), and f^ k \ti) is generated by stochastic recursion 

f ik) (U) = /^(ii-i) + n-WTijrn, (4.1) 

where (r]i) is a white noise, independent of (£»), with Er]\ = 0, Er\\ = 1 and 
7 is an arbitrary nonzero parameter. For the observation model 

*i = / (0) (*i-i)+<T& 
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we apply the estimator given in (p.l|). The resulting errors 5{ = F(ti)—F(ti), 
i = 1, . . . , are defined by a recursion 



1 



5i = 5i-i + -a<5i_i + q„cr£i - q n A5i-\ 
n 



n 2 f+ 1 b7?7j. 



Then, for A, = C n 5i, i > 1 we obtain 

2/3 . . P 

Aj = Aj_i + n 2/3+ 1 (a - q A) A,_i + n 2 ' 3+1 qa^ - n 2 f+ 1 bjrji 
and supply A = 8 . Denote Q? = EA,A*. From (O) it follows 



(4.2) 



■2/3 



I + n 2/5+1 (a - qA) )Q^Al + n 2 /3+i ( a - qA 



2/3 



2/3 9 t 2/3 

+n 2,3+1 (j^qq* + n 2/3+1 ^bb 

2/3 „ 2,3 
= DnQ^Dl + n~^a 2 qq * + n~^-f 2 bb*. 

Similar to (M™, P") (see previous section), let us introduce Q s l : 

+ n"27m ( G - qA)Q^_ 1 (0 - qA)* + cr 2 qq* + 7 2 bb* W. 



Applying the Arzela-Ascoli theorem technique it can be readily shown that 
converges in the local uniform topology to Q s , where 



and lim, 



Qs = Qo + J ((a - qA)Q s , + Q s ,(a - qA)* + a 2 qq* + 7 2 bb*)ds, 

Q s := Q with Q being the unique solution of Lyapunov equation 

(a - qA)Q + Q(a - qA)* + <r 2 qq* + 7 2 bb* = 0. (4.3) 

The matrix Q is a function of arguments q and 7: Q = Q(q,7). We choose 
q = q(7) so that for any 7 

Q(q,7) >Q(q(7),7) ■= Q(l) > 0- (4.4) 
Due to the Kalman filtering theory, the lower bound ( |4.4| ) holds true for 

Q{l)A* 



q(7) 



(7- 



with (3(7) being solution of the algebraic Riccati equation 



aQ(7) + Q(7K+7 2 bb* 



Q{i)A*AQ{i) 



0. 



(4.5) 



(4.6) 



It is well known (see e.g. Theorem 16.2 in Q) that (4.6) possesses a unique 
positive-definite solution provided that block-matrices 

( A \ 

Aa 

Gi = and G 2 =(bb* abb* ... a k bb*) 

\Aa k j 

have full ranks r = k + 1. Notice that G\ is a unite matrix and the rank of 
G2 is k + 1. Consequently, the eigenvalues of the matrix (a — q(j)A) with 
q(7) defined in (fD3) have negative real parts (see, Lemma 16.11 in M). 
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4.3. Minimization of the cost function. The one parameter family 

r={ 7 > 7E >0:,( 7 ) = «^ 



permits a simple numeric implementation and guarantees filtering stability 
mentioned above. In this class we use a constrain parameter / j £ to compen- 
sate unacceptably large boundary layer when the minimization procedure 
yields small values of 7 = argmin 7>0 C(q(7)), with C(q("f)) given in ( 3.11 ). 



The magnitude of 7 e is dictated by k, a, initial conditions and boundary 
layer specifications. The minimization in our class provides 

7 = argminC7(q( 7 )) and q( 7 °) = Q ^ A * . (4.7) 

7>7e a 



5. Explicit minimization procedure 

5.1. Filter gain. In this section we describe a structure of q( 7 °). Recall 
that q( 7 ) = and Q( 7 ) solves the Riccati equation ( |4.6| ) for any fixed 

a. For notational convenience replace Q( 7 ) by Q(7, a) and set U = Q(l, 1). 
Clearly, U solves the algebraic Riccati equation 

aU + Ua* + bb* - UA*AU = 0. 

Kalachev || shows that 

2 l\ fc+l 



a 



Qij(l,<r) = Uijcr -) , i, j = 0, 1, . . . ,k, 



where <5ij(7,cr) and Uij are entries of Q(j, a) and U respectively. Hence, 

l/fc+i 



/j\2/k+l 

91(7) = u oiy-) 



For k < 4, these values are given in the table below. 



k 




U01 


U02 


U 03 


C/04 





1 


NA 


NA 


NA 


NA 


1 


V2 


1 


NA 


NA 


NA 


2 


2 


2 


1 


NA 


NA 


3 


\/4 + \/8 


2 + V2 




1 


NA 


4 


1 + V5 


3 + V5 


3 + V5 


1 + V5 


1 



(Table 1) 



The complex structure of C(q(7)) does not provide an insight of 7 and 
L connection. Numerical simulations show that for a wide range of values 
log(7°) is almost proportional to log(L) (see also Figure |2|). This remark 
enables to construct a simple interpolation tables for the values of 7 and 
C(q(7)) with respect to the parameter L. 
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o = 0.25 



o=1 



10 



10 



~io 4 



10' 



10" 



10"' 



10" 



10 



:io 



10 



10" 



10" 



10" 10' 

7 

o = 4 



10 10 

7 



10' 



10 



10 



10 



~io 4 



10' 



10" 



10" 



10" 



10' 



10' 



L=1 

L=10 

L=100 



L=1; o=0.25; y=0.74082; C=5.278; 
L=10; a =0.25; y=4.4817; C=102.1574; 
L=100; o=0.25; Y=24.5325; 0=2695.7356; 
L=1; o=1; y=1; C=1 8.8975; 
L=10; o=1; y=6.0496; C=260.7145; 
L=100; o=1; y=33.1155; 0=5839.3352; 
L=1; o=4; Y=1.3499; C=92.0461; 
L=10; o=4; y=8.1662; 0=787.7197; 
L=100; o=4; y=49.4024; 0=13850.9423; 



Figure 2. Cost function C(q(7)) in logarithmic scale for 
various L and a; k = 2. 



5.2. Eigenvalues of (a — q(-y)A). Although the eigenvalues of q(7) have 
negative real parts, we may not formally guarantee that they are different. 
So, for k < 4 we give the eigenvalues: 



k = 
k = 1 

k = 2 
k = 3 
fe = 4 



7 s 1/4 



a 

(l; 0.809 ± i0.588; 0.309 ± *0.95l) . 



0.924 ± i0.383; 0.383 ± i0.924 



For k > 4, the fact that the polynomial has different roots should be verified. 

Notice that (7 / 'a) 1 1 is a natural control parameter defining the size 
of boundary layer and should be limited from below by (%/&) 1 '( k+1 ' with 
appropriate j £ . 
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5.3. Cost function. The index form of (a — qA)M(q) = —Lb i.e. 



/ 



-qo 



1 
10 



Qk-i 
\-q k 

allows to find the solution 

1 



0\ 


1 



I Mo(q) \ 
Mi(q) 

Mfc-i(q) 
V M M J 



L 



M 
o 



vv 



M (q) = L-,M 1 (q) = L^,... ,M fc (q) =L 



qk 



qk 



Thus, trace ( M(q)M*(q) J = L 2 ( i ) +£i=o ( f ) . FYom Q it follows 



fc-l ( q 3 



P(q) = a 2 / °° (Ce^-^'e^-^qdt. 

As a result, the final expression for the cost function is 



q* e (a-q^)** e (a~q^)* 



qdt + I-) 2 



^) 2 + £E 



fc-i 



gfe 



3=0 



qj 



qk 



6. Conclusion remark 

In this paper, we use the fact that a class of Kalman filters, being adapted 
to a nonparametric statistic setting, provides the optimal rate of convergence 
in sample size (n — ► oo). We show how to evaluate a normalized risk function 
for large sample size and minimize that value in some subclass of Kalman 
filters with constant filter gain. The Kalman type estimator, as any on- 
line estimator, has inevitable boundary layer. We suggest to reduce the 
boundary layer by interpolation procedure and limitation from below for 
filtering gain. 
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